Chaotic dynamics of a classical radiant cavity 
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The statistical properties of a classical electromagnetic field in interaction with matter are nu- 
merically investigated on a one-dimensional model of a radiant cavity, conservative and with finite 
total energy. Our results suggest a trend towards equipartition of energy, with the relaxation times 
of the normal modes of the cavity increasing with the mode frequency according to a law, the form 
of which depends on the shape of the charge distribution. 
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The study of quantum deviations from classical ergod- 
icity has occupied much of Quantum Chaology since its 
origins. Remarkably enough, the historical development 
of quantum mechanics started with the Blackbody prob- 
lem, which displays a deviation as blatant as possible 
from classical ergodicity. When a classical radiation field 
interacts with matter inside an enclosure with perfectly 
reflecting walls, approach to statistical equilibrium - if at 
all possible - appears to entail unending escape of energy 
towards higher and higher frequencies, in sharp contrast 
to the Planck distribution law Thus the problem 
of blackbody radiation was at once the first problem in 
quantum mechanics, the first problem in quantum field 
theory, and the first problem in quantum chaos. Con- 
siderable progress has been meanwhile attained in un- 
derstanding the complex behaviour of nonlinear classical 
dynamical systems and its quantum counterparts, so a 
re-examination of the blackbody problem in the light of 
such developments appears necessary. One would like, 
first, to build a Hamiltonian model of a radiant cavity, 
which does indeed exhibit the sort of tendency to equipar- 
tition expected in Jeans's time; second, to understand 
how does Planck's law emerge from the quantal dynam- 
ics of that very model. Neither of these issues seems to 
have been satisfactorily dealt with as yet. In this Letter 
we accomplish the first part of the task by presenting a 
model whose classical dynamics leads to equipartition in 
the sense of Jeans. 

Our model is a variant of one which was introduced 
years ago to this purpose, and which was later inves- 
tigated in several papers ||^ . None of those investigations 
was able to detect a tendency towards energy equiparti- 
tion among the normal modes of the cavity. Analogous 
results were obtained on different models such as a 
one-dimensional linear string interacting with nonlinear 
oscillators etc. 

The general picture which emerges from all these 
works, in which the Newton-Maxwell equations were nu- 



merically solved, is that there is no tendency to energy 
equipartition among the field normal modes. The rea- 
son lies with two fully general aspects of the classical 
field-matter interaction. First, the total energy is fi- 
nite, whereas the number of freedoms is infinite; sec- 
ond, field modes can only exchange energy via interaction 
with finitely many mechanical freedoms. Thus mechani- 
cal nonlincarities become less and less effective as energy 
flows from the matter to the field: this fact prevents the 
appearance of an altogether chaotic dynamics, thus caus- 
ing high frequency modes to be nonergodically "frozen" . 
One therefore needs a model, giving rise to chaotic be- 
havior of the mechanical freedoms, no matter how small 
their energy is. 

Though extremely simplified, our model displays this 
property. Let us first consider [|| an electromagnetic 
field confined inbetween two parallel, perfectly reflecting 
plane mirrors, a distance 21 apart. We take Cartesian 
coordinates XYZ with the X axis normal to the mir- 
rors, and restrict to excitations only dependent on X, 
thus getting a 1-dimensional radiant cavity, the normal 
modes of which have angular frequencies LOn = {nc/2l)n, 
n = 1,2,.... Then we introduce a uniformly charged, 
infinite plate of thickness 25, situated midway between 
the mirrors and parallel to them, bound to move along 
the Z-direction only. We denote z its displacement in 
that direction, a and m the charge and mass densities 
per unit surface of the plate, f{x) the normalized (trans- 
verse) distribution of charge in the plate. Finally, the 
plate is subject to a mechanical restoring force per unit 



surface, F{z) 



\z. Using the Coulomb gauge, plus 



zero boundary conditions on the mirrors for the Z com- 
ponent of the vector potential, we obtain the following 
Hamiltonian for the full system plate plus field: 
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where {z,pz) and (qmPn) are canonical conjugated vari- 
ables for the plate and the n— th mode of the field re- 
spectively. In particular, qn{t) is the amplitude of the Z 
component of the vector potential on the n-th normal 
mode of the free field, and the coefficients a„ are given 
by a„ = J_gdxf(x)cos{u}nx/c). means the sum over 
odd n's only, because, with the chosen boundary condi- 
tions, even modes do not interact with the plate. 

Finite-energy states of our hamiltonian system cor- 
respond to vectors in the Hilbert space Tig of square- 
summable, oo-dimensional vectors 



q = 



{ml 
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whose squared norm is just twice the energy. Since 
Hamilton's equations are linear, the evolution of states 
in Tio is given by a unitary group exjp(iHot), with a com- 
plete set of orthonormal eigenvectors , 
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where e = 2cr(7r/mZ)i/^ and Ck is a normalization con- 
stant. These eigenvectors define normal modes of the 
total system, with eigenfrequencies ilu given by the imag- 
inary roots of the secular equation: 



n = l 



0. 



(4) 
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We shall now introduce a nonlinear mechanism, which 
will couple these normal modes, giving rise to energy ex- 
changes between them. To this end we introduce a sec- 
ond plate, parallel to the first, bound to move in the 
same direction, and with the same mass density. We as- 
sume for simplicity this plate to carry no charge, so that 
its motion is not infiuenced by the field. The interac- 
tion between plates is purely mechanical, and simulated 
by elastic bounces when \z{t) — zi{t)\ = R, with zi the 
Z-coordinate of the neutral plate. 

Between collisions, the system is intcgrable: the mo- 
tion of the charged plate and the field is given by the 
action of the group exjp(iHQt), while the second plate 
is moving at a uniform speed. This is most easily de- 
scribed by adding to the vector q one more component 
im}/'^w, with w the velocity of the uncharged plate. The 
new Hilbert space Ti of such q-vectors is the phase space 
of our model. The evolution between collisions is again 



unitary, given by exp{iHt), with a generator H, and a 
complete set of normal modes, which are trivially related 
to the above described ones. At collisions the two plates 
exchange their velocities, thus mixing all the amplitudes 
in the expansion of the state vector over normal modes. 

The evolution from immediately after one collision to 
immediately after the next is given by a map, which, in 
Hilbert space notations, has the following simple form: 



5(q) = {Id-2P)e 



iHt{q) 



(5) 



The 1st (operator) factor describes a collision: Id is 
the identity operator, P is a one-dimensional projection: 
Pq =< e|q > e, where e is the unit vector such that the 
scalar product < e|q > yields the relative velocity of the 
two plates in the state q. The 2nd factor describes evo- 
lution over the free-flight time t{q), which is the smallest 
positive root of the equations 



z{e'"'q) - z(q) - w{q)t\ = 0, 2R. 



(6) 



The map allows for efficient numerical simulation: e.g., 
in the case of 200 oscillators, we were able to follow a 
trajectory up to 5 x 10^ bounces with a relative error in 
energy conservation less than 10^^". 

In numerical simulations, one has of course to consider 
a finite number N of field oscillators. In our computa- 
tions we have varied N and all other parameters, except 
? = 7r,TO = l,c=l; moreover, since the dynamics de- 
pends on energy only via the scaled parameter R/^/E, 
we have always taken E ^ 1 and varied the "free path" 
R instead. 

The choice of the charge density f{x) is important, 
because the coupling of individual modes to the charged 
plate is scaled by the coefficients a„ of the Fourier ex- 
pansion of f{x). Choosing a singular density f{x), 
as in earlier studies, results, at all times, in a power 
law decay of the distribution of energy over the field 
modes, so that truncation effects are already significant 
at small integration times. We have therefore chosen 
f{x) = kexp{—S'^/{6^ — x^)) (the standard compactly 
supported C°° function), with the constant k fixed from 
normalization; this ensures a faster than algebraic, albeit 
nonexponential, decay of the distributions. 

Even though we do not have rigorous results, the dy- 
namics of this billiard-type model appears to be com- 
pletely chaotic independently of the total energy. More- 
over, the finite-dimensional reduced dynamics has posi- 
tive maximal Lyapunov exponents A, Ac (the former being 
defined with respect to real time, the latter to the num- 
ber of collisions). These were numerically computed by 
multiplying matrices obtained from linearization of the 
map (|^) along a trajectory. The exponent A decreases 
with the number of normal modes taken into account, 
because bounces become less frequent, the two plates go- 
ing to rest, in time average, for iV — > oo (see below). On 
the contrary Ac was observed to increase with N; indeed, 
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as collisions become more distant in time, phases change 
more drastically in between them, and their randomiza- 
tion is faster. If instead N is increased keeping the energy 
per mode E/N fixed, both A and appear to saturate, 
suggesting that Lyapunov exponents converge to a finite 
non zero value in the thermodynamic limit. The maximal 
Lyapunov exponents remain positive on reducing cr, with 
no stochasticity threshold displayed. However, the time 
required to reach a converged value becomes larger, be- 
cause trajectories need more time to fill the phase space. 
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FIG. 1. Time-average kinetic energy T of the charged os- 
cillator as a function of the number A'^ of odd field modes 
taken into account. Initial conditions: {z,z-i_) = (0,0), 
— (^1,^2), with ^1,1^2 random numbers such that 
= E = 1, {qn,Pn) = (0,0) Vn = l,2,...,iV; 
coupling constant a — 1, billiard dimension R = 1, frequency 
oJo = 5, thickness parameter 5 = 0.05. Numerical data follow 
the equipartition values T = 1/{2N + 3) (full line). 

Most of our numerical experiments were meant to un- 
derstand how the energy is distributed among all the de- 
grees of freedom (in time average). In Fig|| the time- 
average kinetic energy T = limt^oo T(t) of the charged 
plate is shown as a function of the number N of field 
modes considered. In the above definition, T{t) is the 
time average up to time t of ^mz'^{t), while the limit 
means that the motion has been followed until stabiliza- 
tion of the time-average. Results are in accordance with 
the equipartition theorem: the total energy is equally 
shared between the 2N + 3 relevant canonical variables. 
For N ^ 00 we can extrapolate T = 0, that is, the elec- 
tromagnetic field acts as a friction force on the plate. 

The approach to equilibrium is not uniform, because 
the relaxation time associated with the n-th overall nor- 
mal mode increases with n. To analyze this increase we 
have used the equipartition indicator 



f - 1 

^eff (t) = exp <^-^£;„(t)ln£;„(OL (7) 

where En(t) indicates the normalized time average en- 
ergy (up to the time t) of the n-th normal mode {Ej^+i 
refers to the energy of the neutral plate). The param- 
eter rieff is a measure of the number of modes signifi- 
cantly excited at time t; if only a finite number of modes 
is considered, it also measures the degree of equiparti- 
tion, because rief / = 1 if only one normal mode is ex- 
cited, whereas the maximal value n^ff ~ TV -I- 3/2 is only 
attained in the presence of complete equipartition. As 
far as the numerical simulation is truly representative of 
the infinite-dimensional system, rig/ / appears to increase 
with t slower than any power, but faster than logarith- 
mically (Fig||). 
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FIG. 2. Equipartition parameter Tie// as a function of 
time, for iV = 50 (dashed line), N = 100 (dotted line), 
iV = 200 (full line), with a = 0.3, R = 1, ujo = 0.7, S = 0.05; 
initially half of the energy is given to the discharged plate, half 
to the first overall normal mode. The insert shows rie// as a 
function of the number of bounces (full line) and its random 
phase approximation (dashed line). 

A quantitative description of the relaxation process 
and analytical estimates of expression (^ can be ob- 
tained from cqn.(^ by a random-phase approximation. 
Let us consider an ensemble of trajectories which start at 
time zero from a single normal mode, with randomly dis- 
tributed phases. Denoting -E„(r) the ensemble-averaged 
(normalized) energy on the n— th normal mode after the 
T— th collision, and assuming complete randomness of 
phases at the collision time, from (|^) we get: 

Ekir) = (1 - 4\ek\^)Ek{T - I) + W(r - l)\ek\^ (8) 

where 
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ek=<u„\e>; W{t) =Y,\ek\''Ek{T). (9) 

k 

Recalling the meaning of e, one easily realizes that W{t) 
is proportional to the average kinetic energy of the rel- 
ative motion of the plates. Eqns. (||,^) can be solved 
numerically, to find how Ugf / increases with the number 
of collisions. The result (shown by the dashed line in the 
insert of Fig.(H)) matches quite well with the mmrerical 
solution of the exact equations of motion, confirming the 
validity of the random phase approximation, hence the 
chaotic nature of dynamics. One can also solve (^,^ ana- 
lytically, by implementing a continuous time approxima- 
tion, plus standard Laplace transform techniques. Omit- 
ting details, one finds that the large-r asymptotics of the 
solution is determined by the large-fc asymptotics of the 
coeflicients at- In case of algebraic decay ak ^ 1^1 
dispensing with prefactors which depend on e, / , c and 
estimating the average time delay between the r -f 1- 
th and the r-th collision as t ^ R/^/W{t), one finds 
"'e//(i) ~ With the charge distribution f{x) 

used in our numerical simulations, we cannot give like- 
wise explicit formulas, due to the complicated decay of 
coefficients ak- However, it is possible to prove that Ueff 
increases with t faster than logarithmically, but slower 
than any power of t, as found in Fig.^. Thus the way the 
relaxation time of modes increases with their frequencies 
is determined by the choice of the charge density. 

To verify that the truncated system numerically in- 
vestigated here really represents (up to a certain time) 
the real, infinite-dimensional, model, we plotted in Fig.^ 
neff as a function of the number N of field oscillators 
taken into account. We found that, at any fixed time, 
neff converges on increasing N, its limit value giving the 
number of overall normal modes significantly excited in 
the infinite-dimensional system. As this value increases 
with time, an equilibrium state is never reached. 

In summary, our numerical experiments have exposed 
a chaotic dynamics in the unusual case of an infinite di- 
mensional, conservative system with a finite total energy. 
The dependence of time-averages on initial conditions 
gets lost, and, for any finite-dimensional reduction, the 
system reaches an equilibrium state, with equipartition 
of energy among the degrees of freedom of the field and 
of the matter. 

From our results we infer that, in the real infinite- 
dimensional problem, there is a trend towards equiparti- 
tion, the finite energy of matter being removed to higher 
and higher frequencies of the field. However, this process 
takes place at a nonuniform rate, as relaxation times of 
normal modes increase with their frequency. Therefore, 
as in an old hypothesis of Jeans ||l|, a real equilibrium 
state is never reached. 

Artificial though it may appear, our model is actu- 
ally the simplest, one-dim. model with charged particles 
undergoing elastic collisions inside a reflecting enclosure. 



We believe that more realistic models displaying the same 
basic features will display a similar behaviour. We also 
submit that quantization of this, or similar, field mod- 
els may open a challenging new direction in the field of 
Quantum Chaos, starting from a very old problem. 
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FIG. 3. Equipartition parameter n^ff as a function of the 
number A'^ of odd field modes considered, at t = 10^ (trian- 
gles), t = 10* (squares), t — 10^ (circles), t = 10* (diamonds), 
with initial conditions and parameters chosen as in Fig.| The 
full line indicates equipartition values Me// ~ N + 3/2. 
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